11 Heat equation
In Chapter 3, we derived a general formulation for a balance law. So far, however, we only exploited the balance of mass and momentum (as well as combinations of both) to derive various types of process models. Now, we will add to that the conservation of energy. Our first mathematical model to derive and investigate will be the heat equation to understand the evolving temperature distribution. Here, we will consider an incompressible situation, in which \(\rho\) can be assumed constant.
The heat per unit mass is given by \[c_p T, \quad \text{with} \quad [c_p T] = \tfrac{J}{kg}.\]
Here, \(c_p\) denotes the specific heat capacity (\(\tfrac{J}{kg K}\)) and \(T\) the temperature in \(K\). The heat per unit volume is then given by
\[\rho c_p T, \quad \text{with} \quad [\rho c_p T] = \tfrac{J}{m^3}\]
where \(\rho\) denotes density. Substitution into the general balance law results in
\[ \underbrace{\partial_t \left( \rho c_p T \right) + \nabla \cdot \left( \rho c_p T \mathbf v\right)}_{\text{total derivative of heat per unit volume}} = - \nabla \cdot\underbrace{\mathbf{J}}_{\text{diffusive flux}} + \underbrace{\mathbf{S}}_{\text{heat source terms}} \]
For time being, we will assume heat capacity \(c_p\) as well as density \(\rho\) to be given as constant material parameters. Still, the resulting process model is not closed, as it contains temperature \(T\), velocity \(\mathbf v\), flux \(\mathbf{J}\) and source \(\mathbf{S}\) as unknowns. Note, that velocity \(\mathbf v\) is no longer an unknown of the system, once we solve for the energy balance in conjunction with the previously discussed mass and momentum equation. Additional closure relations are needed!
11.1 Fourier’s law
Newton made the following observation: The rate of heat loss in a one-dimensional bar is proportional to the temperature gradient with a negative sign!
Under the assumption of an isotropic material and assuming that the heat flux in one direction does not affect the heat flux in another direction, Newton’s observation translates into so-called Fourier’s law
\[ \mathbf{J} = - \kappa \nabla T. \]
Here, \(\kappa\) denotes the thermal conductivity measured in \(\tfrac{W}{K m}\). The full energy balance then reads
\[ \partial_t \left( \rho c_p T \right) + \nabla \cdot \left( \rho c_p T \mathbf u\right) = \nabla \cdot \left( \kappa \nabla T \right) + \mathbf{S} \]
In combination with the incompressible Navier-Stokes model this is known as the Navier-Stokes-Fourier system. It models simultaneous incompressible flow and temperature evolution and reads: \[ \begin{aligned} \nabla \cdot \mathbf v &= 0 \\ \partial_t \mathbf v + \mathbf v \cdot \nabla \mathbf v &= - \frac{1}{\rho}\nabla p + \nu \triangle \mathbf v + \mathbf g\\ \partial_t \left( \rho c_p T \right) + \nabla \cdot \left( \rho c_p T \mathbf v\right) &= \nabla \cdot \left( \kappa \nabla T \right) + \mathbf{S} \end{aligned} \]
11.2 Boussinesq approximation
Many relevant situations correspond to a physical regime, in which density variations with temperature cannot be neglected per se, but depend only weakly on a change in temperature (e.g. in water and air moving at moderate velocities). In such regimes, it makes sense to consider density gradients to vanish, while density itself might be subject to a small variation. Very famous is to assume the following linear dependency referred to as the Boussinesq approximation:
\[ \rho = \rho_0 \left( 1 - B (T-T_0) \right) \]
Coefficient \(B\) refers to as the Boussinesq number. Together with the above assumption, we yield the following variant of the momentum equation:
\[ \rho_0 \left( \partial_t \mathbf v + \mathbf v \cdot \nabla \mathbf v \right) = - \nabla p + \mu \triangle \mathbf v + \rho_0 \left( 1 - B (T-T_0) \right) \mathbf g. \]
With that, the full Navier-Stokes-Boussinesq-Fourier system in a gravitation align coordinate system reads
\[ \begin{aligned} \nabla \cdot \mathbf v &= 0 \\ \partial_t \mathbf v + \mathbf v \cdot \nabla \mathbf v &= - \frac{1}{\rho_0}\nabla \left( p - \rho_0 g z \right) + \nu \triangle \mathbf v - \mathbf g B (T-T_0)\\ \partial_t \left( \rho c_p T \right) + \nabla \cdot \left( \rho c_p T \mathbf v\right) &= \nabla \cdot \left( \kappa \nabla T \right) + \mathbf{S}. \end{aligned} \]
Inspection of the system shows, that the temperature couples back into the vertical momentum equation via a source term. This for instance induces well known thermal plumes in the atmosphere.
11.3 Classical heat equation
Let’s consider an idealized situation, in which we have no heat production or decay, \(\rho, \kappa\) and \(c_p\) are constant and \(\mathbf v=0\). This yields what is known as the classical heat equation
\[ \begin{aligned} \rho c_p \partial_t T &= \kappa \triangle T. \end{aligned} \]
The heat equation is often written in terms of the thermal diffusivity \(\alpha = \tfrac{\kappa}{\rho c_p}\) ([\(\alpha\)] = \(\tfrac{m^2}{s}\)):
\[ \begin{aligned} \partial_t T &=& \alpha \triangle T. \end{aligned} \]
A characteristic timescale for diffusive processes is given by \(t_0 = \tfrac{x_0^2}{\alpha}\), in which \(x_0\) stands for the characteristic length scale of interest. For water, we get for instance \(\alpha \approx 10^{-6} \tfrac{m^2}{s}\).
Introducing the scaling
\[ t=t_0 \tilde t = \tfrac{x_0^2}{\alpha} \tilde t \qquad x = x_0 \tilde x \qquad T = T_0 \tilde T \]
allows us to write the heat equation in dimensionless form:
\[ \begin{aligned} \partial_{\tilde t} \tilde T &= \tilde{\triangle} \tilde T \end{aligned} \]
11.4 Fundamental solution
We now would like to derive an analytical solution (that we could potentially test a numerical scheme against). Let’s consider an infinite domain in 1D, namely \(x \in [ -\infty, \infty]\), with \(\lim_{x \to \pm \infty} T = 0\). Note, that energy conservation implies
\[ \int_{-\infty}^{\infty} T dx = const. \]
A solution to the one dimensional heat equation \[ \begin{aligned} \partial_t T &= \partial_x^2 T \end{aligned} \]
can be constructed from a so-called fundamental solution. To understand those, we imagine that heat is released at one specific location \(\zeta\) and at a certain time \(\tau\). This can be represented by the Dirac delta function
\[ \delta(x) = 0 \quad \text{if} \quad x \neq \zeta \quad \text{and} \int_{-\infty}^{\infty} \delta(x) dx = 1. \]
The solution to the heat equation is then given by
\[ G(t,x;\tau,\zeta) = \begin{cases} 0 & t < \tau \\ \delta(x) & t = \tau\\ \frac{1}{2 \sqrt{\pi (t-\tau)}} e^{\frac{- (x-\zeta)^2}{4(t-\tau)}} & t > \tau, \end{cases} \]
which can be verified by substitution into the heat equation. Also, we find \(\lim_{x \to \pm \infty} G(t,x;\tau,\zeta) = 0\), so the far field boundary condition is fullfilled. Self exercise!
This fundamental solution allows us to construct more general solutions via a convolution, and exploiting the superposition principle. This is possible due to the linearity of the process model. For an initial profile \(T_0(x)\) at time \(\tau = 0\), we get
\[ T(t,x) = \int_{-\infty}^{\infty} G(t,x;0,\zeta) T_0(\zeta) d \zeta. \]